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Abstract — We present a method for nonlinear parametric 
optimization based on algebraic geometry. The problem to 
be studied, which arises in optimal control, is to minimize a 
polynomial function with parameters subject to semialgebraic 
constraints. The method uses Grobner bases computation in 
conjunction with the eigenvalue method for solving systems of 
polynomial equations. In this way, certain companion matrices 
are constructed off-line. Then, given the parameter value, an 
on-line algorithm is used to efficiently obtain the optimizer of 
the original optimization problem in real time. 

I. INTRODUCTION 

Optimal control is a very active area of research with 
broad industrial applications [1]. It is among the few control 
methodologies providing a systematic way to perform non- 
linear control synthesis that handles also system constraints. 
To a great extent, it is thanks to this capability of dealing 
with constraints that model predictive control (MPC) has 
proven to be very successful in practice [2], [3]. 

Model predictive control uses optimization on-line to 
obtain the solution of the optimal control problem in real 
time. This method has been proven most effective for 
applications. Typically, the optimal control problem can 
be formulated into a discrete time mathematical program, 
whose solution yields a sequence of control moves. Out of 
these control moves only the first is applied, according to 
the receding horizon control (RHC) scheme. 

The optimal control problem is formulated as a mathe- 
matical program, which can be a linear program (LP), a 
quadratic program (QP) or a general nonlinear program 
(NLP). For hybrid systems, the corresponding mathematical 
programs can be mixed integer programs - MILPs, MIQPs 
or MINLPs [4]. The class of the optimization problem 
depends on the objective function and the class of systems 
one wants to derive an optimal controller for. 

Technology and cost factors, however, make the imple- 
mentation of receding horizon control difficult if not, in 
some cases, impossible. To circumvent these issues, the 
solution of the optimal control problem is computed off- 
line, by solving the corresponding mathematical program 
parametrically [5]. That is, we compute the explicit formula 
giving the solution of the program (control inputs) as a 
function of the problem parameters (measured state). The 
solution then is efficiently implemented on-line as a lookup 
table. 
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In the present work, we extend the concept of the explicit 
solution to the class of nonlinear polynomial systems with 
polynomial cost function. By polynomial systems we mean 
those systems, whose state update equation is given by 
a polynomial vector field. For this class of systems, the 
resulting mathematical program is a nonlinear (polynomial) 
parametric optimization problem. 

While the explicit solution is not generally possible 
in the nonlinear case, we stress the fact that a partial 
precomputation of the optimal control law is still feasible 
using algebraic techniques [6]. In this paper, we use the 
eigenvalue method [7] in conjunction with Grobner bases 
computation to perform nonlinear parametric optimization 
of polynomial functions subject to polynomial constraints. 

II. PARAMETRIC OPTIMIZATION 

Let u G M. m be the decision-variable vector and x G 
W 1 be the parameter vector. The class of optimization 
problems that this paper deals with can generally assume 
the following form: 



min J(it, x) 



s.t. g 



(u, x) < 0, 



(1) 



where J(u, x) G R[xi, . . . , x n , u±, . . . , u m ] is the objective 
function and g S M[x\, ■ ■ . , x n , U\, . . . , u m ] q is a vector 
polynomial function representing the constraints of the 
problem. By parametric optimization, we mean minimizing 
the function J(u, x) with respect to u for any given value 
of the parameter x G X C R n , where X is the set of ad- 
missible parameters. Therefore, the polynomial parametric 
optimization problem is finding a computational procedure 
for evaluating the maps 



where 



u*(x) 



J*(x) 



J* 



(2) 



X I — ► J*, 

argmin J(u, x) 

u 

min J{u, x). 



(3) 



For the sake of simplicity, we assume that the feasible set 
defined by g(u, x) is compact, therefore the minimum is 
attained. Also, in order for (2) not to be point-to-set maps, 
we focus our attention to one (any) optimizer. 



A. Posing the problem 

Our point of departure is the observation that the 
cornerstone of continuous constrained optimization are 
the Karush-Kuhn-Tucker (KKT) conditions. All local and 
global minima for problem (1) (satisfying certain constraint 
qualifications) occur at the so-called "critical points" [8], 
namely the solution set of the following system: 

V u J(u,x) + ^ q i=1 HiV u gi(u,x) = 

Hi9i{u,x) = 

Hi > y ' 

g(u,x) < 0. 

For the class of problems we consider, the two first relations 
of the KKT conditions (4) form a square system of poly- 
nomial equations. Various methods have been proposed in 
the literature for solving systems of polynomial equations, 
both numerical and symbolic [9], [10], [11]. Here we 
consider symbolic methods since our aim is to solve the 
optimization problem parametrically. We should point out 
that the underlying philosophy is that we aim at moving as 
much as possible of the computational burden of solving 
the nonlinear program (1) off-line, leaving an easy task for 
the on-line implementation. 

B. Off-line vs. on-line computations 

The explicit representation of the optimal control law as 
a state feedback has been successfully investigated for the 
linear, quadratic and piecewise affine case. Among other 
advantages of the explicit representation is that one is 
able to analyze the controller, derive Lyapunov functions 
[12], perform dynamic programming iterations [13] in an 
effective way, even compute the infinite horizon solution 
for certain classes of constrained optimal control problems 
[14]. 

Unfortunately, such an explicit representation is not al- 
ways possible. The enabling factor in the case of linear 
systems (or piecewise affine systems) is the fact that the 
KKT system (4) can be solved analytically. In the general 
polynomial case studied here, we have to solve a system of 
(nonlinear) polynomial equations. The next best alternative 
then to an explicit solution is to bring the system in such a 
form, so that once the parameters are specified, the solution 
can be extracted easily and fast. 

III. THE EIGENVALUE METHOD 

In this section we briefly describe the method of eigenval- 
ues ([7], Chapter 2, §4) for solving systems of polynomial 
equations. This method is used in conjunction with Grobner 
bases to perform parametric optimization. 

A. Solving systems of polynomial equations 

Suppose we have a system of m polynomial equations 
fi in m variables Ui 

fi(ui, . . . ,u m ) = 

(5) 

fm(ui,.-.,U m ) = 0. 



These equations form an ideal / £ K[u\, . . . ,u m ], where 
K denotes an arbitrary field: 

/:=</!,...,/„). (6) 

The solution points we are interested in are the points on 
the variety over the algebraic closure K of K, 

V(I) = {s£ K m : h(s) = 0, . . .J m (s) = 0}, (7) 

i.e. the set of common zeros of all polynomials in the ideal 
I. These points can be computed by means of Grobner 
bases. An obvious choice would be a projection-based algo- 
rithm by means of lexicographic Grobner bases, see ([15], 
Chapter 2, §8). Since the computation of a lexicographic 
Grobner basis is very time consuming, we focus on a 
different method. 

The first step we take towards solving (5) is computing a 
Grobner basis with an arbitrary term-order, e.g. graded re- 
verse lexicographic term-order. We define G = {71 , . . . , 74} 
to be this Grobner basis of I. 

B. The generalized companion matrix 

Consider a polynomial function h £ K[u\, . . . , u m }. The 
Grobner basis G and the division algorithm make it possible 
to uniquely write any polynomial h £ K[m, . . . , u m ] in the 
following form: 

Q 

h = c 1 (u)-f 1 -\ \-c t (u)j t + h , (8) 

where h G is the unique remainder of the division of h with 
respect to the Grobner basis G. The polynomial h can in 
turn be multiplied with another polynomial function / £ 
K[u\, . . . , u m ] and their product expressed as follows: 

/ • h = di(«) 7 i + • • • + dt(u)7t +T~h G ■ (9) 

In the generic case, the ideal I will be zero-dimensional, 
which means that the corresponding quotient ring 

A = K[ Ul ,...,u m ]/I (10) 

is a finite-dimensional X-vector space ([15], Chapter 5, 
§2). The quotient ring of an ideal can be thought of as the 
set of all polynomials that do not belong to the ideal but 
belong to the underlying ring. Denote with b = [bi, ... , 6;] T 
the vector of the standard monomials. A monomial is 
standard if it is not divisible by any leading monomial of a 
polynomial in the Grobner basis. These standard monomials 
of G form a basis 

B = {6i,...,6j} (ID 

for the if -vector space A. As a result, every remainder can 
be expressed with respect to this basis as an inner product 

rt = af-b, (12) 

where dj £ K l . We can now define the map nih : A — > A 
as follows: if p G £ A, then 

Q 

m h {p G ) :=7Tp G = h G -p G £ A. (13) 



The following proposition holds. 

Proposition I: Let h G K[u\, . . . ,u m ]. Then the map 
nih ■ A — > A is if -linear. 

The proof of Proposition 1 can be found in ([15], p. 51). 
Since A is a finite-dimensional vector space and the map 
rrih is linear, its representation with respect to a basis of 
this vector space is given by a square matrix Mh. The I x l- 
matrix Mh is called the generalized companion matrix. 

C. Computing the companion matrix 

To compute the matrix Mh, assume that we have the basis 
B = {bi, . . . , b{\ consisting of the standard monomials 6j of 
the Grobner basis G. Then, for each one of them, compute 
the remainder ri of the polynomial h ■ bi with respect to the 
Grobner basis G: 

-G 



h ■ bi 



Vb t eB. 



(14) 



All ri e A can in turn be expressed as an inner product 

n = af ■ b (15) 

with respect to the basis B. By collecting all vectors a» for 
all basis elements [7], we can construct a representation of 
the map m h with respect to basis B, i.e. calculate the matrix 
M h as follows: 



M h = [oij] 



(16) 



Computing the companion matrix is a standard alge- 
braic procedure implemented in various packages, e.g. in 
Maple 10. 

D. Evaluating polynomial functions on a variety 

Consider a polynomial function h G M\u\, . . . , u m \. The 
amazing fact about the matrix Mh is that the set of its 
eigenvalues is exactly the value of h over the variety V(J) 
defined by the ideal I. More precisely, V(I) is the set of 
all solution points in complex m-space C m of the system 
(5). The following theorem holds. 

Theorem 1: Let I G C[iti, . . . ,u m ] be a zero- 
dimensional ideal, let h G C[ui, . . . , u m \. Then, for A G C, 
the following are equivalent: 

1) A is an eigenvalue of the matrix Mh 

2) A is a value of the function h on the variety V(/). 
The proof can be found in ([7], p. 54). 

To obtain the coordinates of the solution set of (5), we 
evaluate the functions 



hi : u 
h m : u 



Ml 



(17) 



on the variety V(I) defined by the ideal I, where u above 
denotes the vector (m, . . . , u m ). This can be done by means 
of the associated companion matrices of the functions hi. 
The following theorem taken from ([9] p. 22) is the basis 
for the calculation of these point coordinates. 



Theorem 2: The complex zeros of the ideal I are the 
vectors of joint eigenvalues of the companion matrices 
M Ul . . .M Um , that is, 

V(I)={(u 1 ,...,u m )GM. m : 

3«£l m Vi : M Ui v = Ui v} 

It has to be noted that any vector-valued polynomial 
function h : W n — ► M can be evaluated over a zero- 
dimensional variety in the same way. 

IV. THE ALGORITHM 

In this section, we present the proposed algorithm, which 
consists of two parts: the off-line part, where the general- 
ized companion matrices for the optimization problem are 
constructed, and the on-line part where this precomputed 
information is used and given the value of the parameter x, 
the optimal solution is efficiently extracted. 

A. Idea 

Under certain regularity conditions, if J* (defined in (3)) 
exists and occurs at an optimizer u* , the KKT system (4) 
holds at u*. Consequently, J* is the minimum of J(u,x) 
over the semialgebraic set defined by the KKT equations 
and inequalities (4). These conditions can be separated in 
a set of inequalities and a square system of polynomial 
equations. The method of eigenvalues for solving systems of 
polynomial equations as described in section III can be used 
for the latter. This method assumes that the ideal generated 
by the KKT system (4) is zero-dimensional. 

By ignoring the inequalities, a superset of all critical 
points is computed and in a second step, all infeasible points 
are removed. Finally, among the feasible candidate points 
those with the smallest cost function value have to be found 
via discrete optimization. By discrete optimization we mean 
choosing among a finite set that point, which yields the 
smallest objective function value. 

B. Off-line Part 

In K[u\, u m , fj,i, fx q ], where K is the field of 
rational functions R(x\ . . . , x n ) in the parameter x, we 
define the KKT ideal 

q 

Ikkt = (V u J(u,x) + y]fiiV u 9i(u,x), Higi{u,x)) (18) 

i=l 

containing all the equations within the KKT-system (4). All 
critical points for the optimization problem (4) and fixed x 
are the subset of real points on the KKT-variety 



Vkkt ^= Vkkt — V(Ikkt) ■ 



(19) 



Using the method described in section III we can compute 
these by means of the generalized companion matrices. 

The algebraic part of the algorithm, i.e the computation 
of the companion matrices can be done parametrically. 
For one thing, one could use Grobner bases computation 
for the ideal Ikkt and try to compute the corresponding 



companion matrices M Ui and M fli directly. Owing to 
the structure of the polynomial equations of the KKT- 
system (18), this problem is very poorly conditioned. The 
difficulties stem from the fact that the ideal Ikkt is by 
construction decomposable. It contains terms like fngi{u, x) 
which lead to a reducible variety V(Ikkt)- 

To overcome this obstacle, we factorize the generators of 
the Grobner basis (i.e. the polynomials appearing in relation 
(18)) and express the ideal Ikkt as an intersection of 
super-ideals Ij.KKT- The super-ideal Ij.KKT denotes the 
ideal constructed by fixing a subset of p active constraints 
gi(u,x) among the set of all q constraints gi(u,x) - see 
(18). The corresponding Lagrange multipliers are denoted 
with jSj. This leads to 

Ij,KKT = ( V u J(u,x) +J2f =1 ^ugt(u,x), 

gi(u,x) } 

with the feasibility inequalities 



> 
< 








(21) 



9i(x,u) 

Therefore, the ideal Ikkt can be expressed as an intersec- 
tion of 9 := $({gi(u,x)}'-_ 1 ) — 2 q super-ideals, where 9 is 
the cardinality of the power set of all q constraints. Namely, 



Ikkt — Ij.KKT 



(22) 



Relations (20) and (21) lead to a large number of super- 
ideals which are much better numerically conditioned than 
the original problem, even though they are not necessarily 
radical. Since many of the sub-varieties V(Ikkt) are 
empty, a Grobner basis computation for each ideal I^kkt 
identifies these infeasible cases in advance and reduces the 
subsequent companion matrix computations tremendously 
by discarding them. 

The number of solutions over K in the non-empty 
sub-varieties Vj.KKT = V(Ij,KKT) can be calculated by 
means of the Hilbert polynomial ([15], Chapter 9, §3). For 
zero-dimensional varieties this polynomial reduces to an 
integer, which is equal to the number of solutions counting 
multiplicity. 

If the sub-variety has only a single solution, the co- 
ordinates Ui of the candidate solution can be computed 
analytically as a rational function of the parameters x. In 
this case, the polynomials in the Grobner basis from a set 
of linear equations in the decision variables that can be 
solved analytically. For all sub-varieties with more than one 
solution, a companion matrix has to be computed. The result 
are companion matrices whose entries are rational functions 
of the parameter x. 

Specialization of the parameters gives a map from the 
field K to the field R of real numbers. If the real parameters 
are chosen generically enough, then the given Grobner basis 
remains a Grobner basis, but for special choices of the 
parameters some trouble may arise. For instance, it may 
happen that a specialization leads to zero denominators. 



To handle this case, comprehensive Grobner bases can be 
used [16]. The parametric computation is guaranteed to be 
correct only if the sequence of leading coefficients of the 
result and the sequence of greatest common denominators 
removed in the computations are nonzero [16]. If ordinary 
methods such as Buchberger's algorithm are used to com- 
pute Grobner bases, these issues have to be kept in mind. 

A summary of the off-line algorithm appears in Algo- 
rithm 1. 

Algorithm 1 Off-line Part: 

Input: Objective function J(x, u) and constraints 

9i{x,u) < 0. 

Output: Set of feasible sub-varieties Vj.KKT with their 
generalized companion matrices M Uji and Mj^, or 
an explicit function u* i for their candidate optimizer. 
1: for all combination of active and inactive constraints 
do 

construct Ij.KKT 

calc. Grobner basis Gj for Ij.KKT 
if Gj =< 1 > then 

discard the super-ideal 
else 

calculate number of solutions of Vj.KKT by means 
of the Hilbert polynomial 
if ^Vj.KKT = 1 then 

Express all u*^ as rational functions in the 
parameter x 
else 

Compute generalized companion matrices M^ Ui 
and Mj.^ for all decision variables Ui 



10: 
11: 

12: 
13 
14: 
15 
16: 



end if 
end if 
end for 



return: Mj )Ui and Mj^ t , resp. u* } l and p,j ti 



C. On-line Part 

In order to evaluate the point coordinates of the KKT 
sub-varieties, we need to compute eigenvectors and eigen- 
values for the companion matrices. Generally, eigenvalue 
computation cannot be done parametrically. The parameter 
x has to be fixed to a numerical value and this computation 
is done on-line. 

Given the precomputed generalized companion matrices 
Mj iUi and Mj^ (resp. an explicit expression for all sub- 
varieties with linear Grobner basis) for all possible feasible 
combinations of active and inactive constraints, the on-line 
algorithm takes the value of the parameters x to compute 
the optimum J* and the optimizer u*. The three main steps 
of the algorithm are: 

1) calculate all critical points 

2) remove infeasible solutions 

3) find the feasible solution u* with the smallest objec- 
tive function value J* = J(u*). 



Since all companion matrices have been computed para- 
metrically, the remaining part that has to be done is linear 
algebra. For every non-empty sub-variety Vj.KKT, a set 
of right eigenvectors {v} is computed for the companion 
matrices Mj. Ui of the j-th sub-variety, see Theorem 2. 
Because all companion matrices for a sub-variety Vj^KKT 
commute pairwise, they form a commutative sub-algebra 
within the non-commutative algebra oflxl matrices, where 
I is the companion matrix dimension (11), see also [7]. 
Therefore, it suffices to calculate the eigenvectors for a 
single arbitrary matrix in this sub-algebra, because they 
all share the same eigenvectors. To avoid computational 
problems, we choose a matrix Mj^ ran d in this sub-algebra 
as a random linear combination of the companion matrices 
associated with the decision variables Mj }Ui , i.e. 

Mj^and = C\Mj jUl -\ VC m Mj. Um + 

+ c m+1 M hfll H + Cm+pMj^ , 

where Cj G M. are randomly chosen. This ensures, with a low 
probability of failure, that the corresponding eigenvalues 
will all have algebraic multiplicity of one ([7], Chapter 2, 
§4). 

The sets of eigenvectors {v}j can now be used to 
compute all candidate critical points and their Lagrange 
multipliers jlj^ for the sub-variety Vj.KKT- To avoid un- 
necessary computations, we first calculate the candidate 
Lagrange multipliers fij.i for each sub-variety Vj,kkt- 
In this way, complex or infeasible candidate points with 
fij.i < for some i can be immediately discarded before 
the candidate optimizers u* i are computed. For all sub- 
varieties with cardinality one, the problem of computing the 
critical points reduces to an evaluation of the precomputed 
functions. 

For all non-discarded candidate solutions, it remains to 
be checked whether they are feasible, i.e. g(u*i,Xi) < 0. 
To achieve that, a set of feasible local candidate optimizers 
S = {u* i } is initially calculated by collecting all fea- 
sible candidate optimizers. After computing the objective 
function value J(v,j A ,x) for all candidate optimizers, the 
optimal solution 

J* = min J(u*,-, x) 

u . 4 es 

and the optimizer 

u* = arg min J(u^ i; x) 

for the optimization problem (1) can be easily obtained via 
discrete optimization over the finite set S. 

A summary of the on-line algorithm can be seen in 
algorithm 2. 

V. OPTIMAL CONTROL APPLICATION 

In this section we fist give a description of the model pre- 
dictive control optimization problem to show the connection 
of parametric optimization and optimal control. 



Algorithm 2 On-line Part: Companion matrices M Ui and 
Mjx i for all non-empty sub-varieties Vj.KKT, resp. explicit 
expression for cardinality one sub-varieties has to be pro- 

vided. 

Input: Value of the parameter x (state measurement taken 

in real time). 
Output: Optimal cost J* and optimizer u*. 
1: for all feasible sub-varieties Vj.KKT with §Vj t KKT > 1 

do 

2: specialize parameter x in M Ui and 

3: calc. a set of common eigenvectors {v} for the 
companion matrix Mj^ ran d 

4: solve Mj.^v — fij A v to obtain the joint-eigenvalues, 
i.e. candidates for fij.i 

5: discard all eigenvectors with corresp. fij.i < 

6: use the remaining eigenvectors to calc. joint- 
eigenvalues of Mj iUi to obtain candidates for u* A 

7: end for 

8: for all feasible sub-varieties Vj.KKT with $Vj t KKT = 1 
do 

9: evaluate fij.i(x) for all i 
10: if 3i : jij,i(x) < then 
11: discard sub-variety Vj^KKT 
12: else 

13: evaluate Uj ^x) 
14: end if 
15: end for 

16: for all evaluated candidate points {u* t }j do 

17: if g k (u* A ,x) > then 

18: discard candidate point u*^ 

19: else 

20: evaluate J(u* ^x) 
21: end if 
22: end for 

23: compare J(u* A ,x) for the calculated candidates u* A 

and choose optimal J* and corresponding u* 
24: return: optimal cost J* and optimizer u* 



A. Nonlinear model predictive control 

Consider the nonlinear discrete-time system with state 
vector x G R™ and input vector u G M m 

x(k + l) = f(x(k),u(k)) (24) 

subject to the inequality constraints 

g(u(k),x(k)) < 0, k = 0, . . . , N , (25) 

where N is the prediction horizon and g G 
R[xi, . . . , x n ,ui, . . . , u m ] q is a vector polynomial function 
representing the constraints of the problem. We consider 
the problem of regulating system (24) to the origin. For 
that purpose, we define the following cost function 

JV-l 

J(U?-\ x ) = ]T L k (x(k),u(k)) + L N (x(N), u(N)) , 

fe=0 



where fTpf -1 := [u(0), . . . , u(N - 1)] is the optimiza- 
tion vector consisting of all the control inputs for k = 
0, . . . , N — 1 and x(0) = x is the initial state of the 
system. Therefore, computing the control input is equivalent 
to solving the following nonlinear constrained optimization 
program 

min J(Uq^ 1 ,x ) 

u 

x(k + l) = f(x(k),u(k)) (26) 
s.t. g(u(k),x(k)) < 0, k = 0,...,N. 

Forming a vector u of decision variables with Uk = u(k) 
and renaming x(0), problem (26) is written in the more 
compact form 



min J(u,x) s.t. g(u,x)<0, 



(27) 



where J(u, x) is a polynomial function in u and x, u G R m 
is the decision variable vector and the initial state x — 
x{Q) € R™ is the parameter vector. This is exactly problem 
(1), a nonlinear parametric optimization problem. Our goal 
is to obtain the vector of control moves u. 

B. Illustrative example 

In this section we illustrate the application of the pro- 
posed method by means of a simple example. The off- 
line algorithm including the algebraic methods and the 
case enumeration (22) have been implemented in Maple. 
A Maple-generated input file is used to initialize Matlab, in 
order to compute the optimizer on-line. 

Consider the Duffing oscillator [17], a nonlinear oscillator 
of second order. An equation describing it in continuous 
time is 



y(t) + 2(y(t) + y(t) + y(t) 3 = u(t), 



(28) 



where y £ R is the continuous state variable and u G R the 
control input. The parameter ( is the damping coefficient 
and is known (here ( — 0.3). The control objective is to 
regulate the state to the origin. To derive the discrete time 
model, forward difference approximation is used (with a 
sampling period of h — 0.05 time units). The resulting state 
space model with a discrete state vector x € R 2 and input 
u e R is 

xi(k + 1) ] 1 h \ xi(k) 

x 2 (k + l) \ ~ [ -h (I -2(h) \ [ x 2 {k) 





+ 



h 



u{k) + 





-hx\(k) 



An optimal control problem with prediction horizon N 
weight matrices 



3, 



Q = 



R 



1 


1 

To 



and state-constraints 



leads to the following optimization problem: 

J* = 

+ Ei= u ( k + i)R<k + i) 



mm 

u(fe),u(fe+l),«(fe+2) 



xi(k+i) 
X2(k+i) 



s.t. \\x(k +j)\\oo < 5 Vj = 1...N. 

Of these twelve constraints there are ten constraints in- 
volving u(k + i), which have to be considered during the 
optimization. As described in section IV the KKT-variety 
will be split in 2^ 9 ^ = 2 10 = 1024 sub-varieties. For all 
of them a Grobner basis needs to be computed. It turns out 
that only 29 of these are feasible, i.e. having a Grobner 
basis different from unity. Only these cases have to be 
further considered in the online algorithm. Among them 
there are 24 sub-varieties Vj t KKT with a linear Grobner 
basis. For these, a closed form expression for the candidate 
optimizers u* A can be computed. For the remaining five 
cases companion matrices have to be computed, requiring 
eigenvalue computation in the on-line algorithm. These sub- 
varieties Vj.KKT have five solutions counting multiplicities, 
i.e. the companion matrices are 5 x 5 matrices. 

The trajectory of the controlled system starting from an 
initial state of xi(0) = 2.5 and x 2 (0) = 1 is shown in 
Figure 1. Figure 2 shows the state-space evolution of the 
controlled Duffing oscillator and its free response without 
the controller. In the uncontrolled case, a weak dynamic 
behavior and a violation of the constraint ^(t) > —5 can 
be observed. 



- controlled 
free response 



- 
-1 




||z(fc+j)Hoo<5 Vj = l...JV 



urn 



Fig. 1. State-space diagram of the Duffing oscillator 



The precomputation of companion matrices and the so- 
lutions i took less than one minute on a Intel Pentium 
3 GHz with 1 GB RAM. The online algorithm needed less 
than 3.5 s to obtain the global optimum even with a naive 
brute-force on-line search algorithm for the minimization 
over the finite set of candidate points. It has to be noted 
that most of the time of these 3.5 s is consumed by the 
evaluation of expressions with the Matlab Symbolic Math 
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Fig. 2. State and input evolution of the controlled Duffing oscillator 

Toolbox. An efficient implementation, in C for instance, 
would be orders of magnitude faster. 

VI. CONCLUSIONS AND OUTLOOK 

The main contribution of this paper is a new algorithm 
for nonlinear parametric optimization of polynomial func- 
tions subject to polynomial constraints. The algorithm uses 
Grobner bases and the eigenvalue method for solving sys- 
tems of polynomial equations, to evaluate the map from the 
space of parameters to the corresponding optimal value and 
optimizer. The algorithm is very general, computationally 
robust and can be applied to a wide range of problems. 

The punchline of the proposed approach is the precompu- 
tation of the generalized companion matrices, thus partially 
presolving the optimization problem and moving the com- 
putational burden off-line. The method has been developed 
with model predictive control in mind. The connection to 
optimal control problems has been illustrated by applying 
the method to the Duffing oscillator. 

Finally, there is ongoing research on exploiting the 
structure of specific control problems, including sparseness 
and genericity assumption relaxation. More specifically, 
sparse resultant techniques are investigated to compute the 
companion matrices. Combining this method with recently 
proposed "Sum of Squares Programming" methods, based 
on semi-definite representations of finite varieties [18], 
seems to be a promising direction for further research. 
Moreover, the integration of the proposed scheme with 
dynamic programming is also explored. 
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